OpenEO Platform - User Story 1

Access of pre-computed data over a large area.

from eo_utils import *

Background - Explore the openEO Platform

Before loading and processing the data, we can have a look at the available collections:

backend = "https://openeo.cloud"
con = openeo.connect(backend)
con.list_collections()

Get the metadata for Sentinel-2 ARD collection:

con.describe_collection("boa_sentinel_2")

Authenticate using EGI

conn = openeo.connect(backend).authenticate_oidc(provider_id='egi')

Area of Interest selection

Run the next cell, click on the rectangle button and then draw the area over the map like in the following example image

center = [47.29, 10.75]
zoom = 7

eoMap = openeoMap(center,zoom)
addS2Tiles(eoMap)
eoMap.map

Get the bounding box from the previous map.

bbox = eoMap.getBbox()
print('west',bbox[0],'\neast',bbox[2],'\nsouth',bbox[1],'\nnorth',bbox[3])

west 11.066666 
east 11.171036 
south 46.038835 
north 46.110286

Sentinel-2 L2A ARD (Analysis Ready Data) access

1. Define collection, time range and bands of interest.

collection      = 'TERRASCOPE_S2_TOC_V2'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2020-08-08", "2020-08-12"]
bands           = ["B04"]

s2_l2a = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)

1. Save operation as final step to the pipeline.

s2_l2a_geotiff = s2_l2a.reduce_dimension(reducer="mean", dimension="t").save_result(format="GTIFF")

1. Create a new batch job on the back-end and ask to process it.

job = conn.create_job(s2_l2a_geotiff)
job_id = job.job_id
if job_id:
    print("Batch job created with id: ",job_id)
    job.start_job()
else:
    print("Error! Job ID is None")

Batch job created with id:  vito-66cc8836-1370-4682-8332-10a8ce9bbc21

We can get a description of the job and check its status.

job = conn.job(job_id)
job_description = job.describe_job()
print("Batch job with id: ",job_id, ' is ',job_description['status'])

Batch job with id:  vito-66cc8836-1370-4682-8332-10a8ce9bbc21  is  running

job

results = job.get_results()
results

1. Once the job status is 'finished', we can download the result and visualize it.

results.download_file("./data/US1/s2_result.tif")

Running the next cell, the Sentinel-2 data will appear on the map as in the following screenshot.

addLayer(eoMap,'./data/US1/s2_result.tif','S2 B04',clip=[0,1500])
eoMap.map

B04_S2 = xr.open_rasterio('./data/US1/s2_result.tif')

B04_S2

<xarray.DataArray (band: 1, y: 829, x: 965)>
[799985 values with dtype=uint16]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 5.156e+06 5.156e+06 5.156e+06 ... 5.148e+06 5.148e+06
  * x        (x) float64 6.743e+05 6.743e+05 6.743e+05 ... 6.839e+05 6.839e+05
Attributes:
    transform:      (10.0, 0.0, 674250.0, 0.0, -10.0, 5156280.0)
    crs:            +init=epsg:32632
    res:            (10.0, 10.0)
    is_tiled:       1
    nodatavals:     (0.0,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

fig, ax1 = plt.subplots(1,1,figsize=(20,20))

ax1.imshow(B04_S2[0].values,cmap='Greys_r',vmin=0,vmax=2000)

<matplotlib.image.AxesImage at 0x7fb4911ee160>

Sentinel-1 ARD access

We use the same area as before following the same steps with a different collection

collection      = 'TERRASCOPE_S1_GAMMA0_V1'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2019-09-18", "2019-09-25"]
bands           = ["VH"]

s1 = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)
s1_geotiff = s1.reduce_dimension(reducer="max", dimension="t").save_result(format="GTIFF")

1. Create a new batch job on the back-end and ask to process it.

s1_job = conn.create_job(s1_geotiff)
s1_job_id = s1_job.job_id
if s1_job_id:
    print("Batch job created with id: ",s1_job_id)
    s1_job.start_job()
else:
    print("Error! Job ID is None")

Batch job created with id:  vito-9879f280-2f9d-40b1-950c-81a5263296a5

s1_job = conn.job(s1_job_id)
job_description = s1_job.describe_job()
print("Batch job with id: ",s1_job_id, ' is ',job_description['status'])

Batch job with id:  vito-9879f280-2f9d-40b1-950c-81a5263296a5  is  running

results = s1_job.get_results()
metadata = results.get_metadata()
results.download_file("./data/US1/s1_result.tif")

WindowsPath('data/US1/s1_result.tif')

Running the next cell, the Sentinel-1 data will appear on the map as in the following screenshot.

addLayer(eoMap,'./data/US1/s1_result.tif','VH',clip=[0,0.1])
eoMap.map

Landsat-8 ARD access

Accessing data from overlapping areas in the original Landsat tiling system.

We can have a look at the original spatial coverage of the ARD Landsat-8 collection, with its four overlapping areas.

We are going to demonstrate that, with the Equi7 grid used in this ARD collection, the data is being reprojected and therefore there are no problems of duplicated data in overlapping areas.

conn_2 = openeo.connect("https://openeo.eodc.eu/v1.0")

itemIter = conn_2.collection_items(
    "boa_landsat_8",
    temporal_extent = ["2019-09-20", "2019-09-22"],
    spatial_extent = [5.0,45.0,20.0,50.0]
)

next(itemIter)

center = [45.85, 11.0]
zoom = 11

eoMap2 = openeoMap(center,zoom)
eoMap2.map

bbox = eoMap2.getBbox()
print('west',bbox[0],'\neast',bbox[2],'\nsouth',bbox[1],'\nnorth',bbox[3])

west 10.809517 
east 10.91526 
south 45.850916 
north 45.9121

collection      = 'boa_landsat_8'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2019-09-20", "2019-09-22"]
bands           = ["B02","B03","B04"]

l8 = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)
l8_nc = l8.save_result(format="netCDF")

1. Create a new batch job on the back-end and ask to process it.

l8_job = conn.create_job(l8_nc)
l8_job_id = l8_job.job_id
if l8_job_id:
    print("Batch job created with id: ",l8_job_id)
    l8_job.start_job()
else:
    print("Error! Job ID is None")

l8_job = conn.job(l8_job_id)
job_description = l8_job.describe_job()
print("Batch job with id: ",l8_job_id, ' is ',job_description['status'])

1. Once the job status is 'finished', we can download it.

results = l8_job.get_results()
metadata = results.get_metadata()
results.download_file("./data/US1/l8_result2.nc")

PosixPath('data/US1/l8_result2.nc')

dc = xr.open_dataset('./data/US1/l8_result2.nc')
dc

<xarray.Dataset>
Dimensions:      (y: 274, x: 310)
Coordinates:
    time         datetime64[ns] 2019-09-21T12:00:00
    spatial_ref  int32 0
  * y            (y) float64 1.424e+06 1.424e+06 ... 1.416e+06 1.416e+06
  * x            (x) float64 4.815e+06 4.815e+06 ... 4.825e+06 4.825e+06
Data variables:
    B02          (y, x) float64 190.0 192.0 181.0 193.0 ... 263.0 272.0 264.0
    B03          (y, x) float64 369.0 359.0 330.0 373.0 ... 548.0 550.0 541.0
    B04          (y, x) float64 226.0 230.0 213.0 235.0 ... 296.0 311.0 296.0
Attributes:
    crs:                PROJCS["Azimuthal_Equidistant",GEOGCS["WGS 84",DATUM[...
    datetime_from_dim:  2019-09-21T12:00:00.000000000

We can visualize the RGB composite as a simple image using Matplotlib. Notice the seamless integration of the overlapping Landsat tiles.

fig, ax1 = plt.subplots(1,1,figsize=(10,10))

ax1.imshow(tone_mapping(dc.B04,dc.B03,dc.B02),cmap='brg')

<matplotlib.image.AxesImage at 0x7fb492481a90>